| title: “broccoli-population-structure” |
| author: “Sean O’Neill” |
| date: “8/2/2021” |
| output: html_document |
| html_document: |
| highlight: pygments |
| theme: united |
| toc: yes |
| toc_float: |
| collapsed: yes |
| smooth_scroll: yes |
| pdf_document: |
| toc: yes |
rm(list=ls())
#get all the names
taxalabels <- read.delim("taxa/taxaBNUM_NAME.txt", header=F, stringsAsFactors=F, row.names = NULL)$V1
#make duplicate names unique
taxalabels <- make.names(taxalabels, unique = T)
#make some colors. http://vrl.cs.brown.edu/color
clist <- list(
sean=c("#fd00fe", "#62ce75", "#473c85", "#24a5f7", "#9a1073", "#4ad9e1", "#284e37", "#cad7d4", "#723521", "#fea53b", "#c00018", "#f7b5be", "#788c3b", "#5920af", "#eb74c3", "#6975fe", "#238910", "#528e8c", "#65f112", "#a57784", "#b3e61c", "#ae56ff", "#d6c951", "#ff7074"),
shawn=c("#fd00fe", "#84bc04", "#1f3e9e", "#e775cc", "#46ebdc", "#256b33", "#82d1f4", "#7c28a4", "#2af385", "#214a65", "#aaa4e1", "#a1085c", "#148fae", "#cadba5", "#943112", "#fa8d80", "#604020", "#f1d438", "#eb1138", "#f7931e", "#270fe2", "#a07d62")
)
# add length of palettes
lengths <- sapply(clist,length)
names(clist) <- paste0(names(clist),"_",lengths)
# par(mar=c(0.2,6,0.2,0))
# par(mfrow=c(length(clist),1))
#
# #for some reason, barplot() doesn't work in .Rmd
# for(i in 1:length(clist))
# {
# {barplot(rep(1,max(lengths)),col=c(clist[[i]],rep("white",max(lengths)-length(clist[[i]]))),axes=F,border=F)
# text(x=-0.1,y=0.5,adj=1,label=names(clist)[i],xpd=T,cex=1.2)}
# }
Phylogenetic Tree
#source('scripts/rename-taxa.sh --infile="BOL21afZAXA00/RAxML/outraxml.raxml.supportFBP" --outfile="BOL21afZAXA00/RAxML/raxml.tree" --repfile="replaceNEWtoOLD.csv"')
#setwd("BOL21afZAXA00/")
out.group <- c("W001_W001_WT_A","W006_W006_WT","W005_W005_WT","W004_W004_WT","W002_W002_WT","W003_W003_WT")
tree <- read.tree(file="RAxML/raxml.tree") -> uniqtree
tree$tip.label <- make.names(tree$tip.label, unique = T)
tree <- root(tree, outgroup = out.group)
tree <- ladderize(tree)
temp <-plot(tree,node.pos=1, font=1, align.tip.label=TRUE ,plot=FALSE)
is_tip <- tree$edge[,2] <= length(tree$tip.label)
ordered_tips <- tree$edge[is_tip, 2]
write.csv(tree$tip.label[rev(ordered_tips)], file="order.csv", row.names=FALSE)
# # Tree for popstructure
# svg("image/raxml.tree.svg", height = 35, width=20, pointsize = 10)
# par(oma=c(0,0,0,10))
# plot(tree,node.pos=1, font=1, align.tip.label=TRUE, cex=1.5, edge.width=4, no.margin=TRUE, label.offset = 0.01 )
# dev.off()
tree.file <- "../raxml.tree_delete.png"
png(tree.file, units = "in", height = 25, width=14, pointsize = 10, res = 150)
# tree.file <- "images/raxml.tree_delete.svg"
# svg(tree.file, height = 25, width=14, pointsize = 10)
par(oma=c(0,0,0,10))
plot(tree,node.pos=1, font=1,show.tip.label = TRUE, align.tip.label=TRUE, cex=1, edge.width=2, no.margin=TRUE, label.offset = 0.01 )
tree.labels <- tree$node.label
tree.labels[as.integer(tree.labels) > 0] <- ""
#node.indices <- which(as.integer(tree.labels)<50)
#tree.labels <- as.character(tree.labels)
nodelabels(tree.labels, frame = "none", adj=c(-.125),cex=.75)
dev.off()
knitr::include_graphics(tree.file)
PCA
#First, manually create a file called order.groups.csv. BOL2137SMinCn01 has an example
#choose a color list from clist.
gclist <- clist$sean_24
pca <- read.delim("PCA/PCA1.txt",skip = 2)
# pca2 <-read.delim("PCA/PCA2.txt",skip = 2)
# pca3 <-read.delim("PCA/PCA3.txt",skip = 0)
# hist(pca3$Eigenvector4)
pca$Taxa <- taxalabels
order <- read.delim("order.groups.csv")
#rearrange pca rows to be like order
pca <- pca[match(order$taxa,pca$Taxa),]
pca$group <- order$n00
pca$color <- gclist[order$n00]
#now make the group color list
gc <- order$color
axis = list(showline=FALSE,
zeroline=FALSE,
gridcolor=rgb(1,1,1), #'#ffff',
ticklen=0,
automargin = TRUE,
titlefont=list(size=13))
fig11 <- pca %>%
plot_ly(
height=1440,
width=1440
)
fig11 <- fig11 %>%
add_trace(
type = 'splom',
dimensions = list(
list(label='PC1', values=~PC1),
list(label='PC2', values=~PC2),
list(label='PC3', values=~PC3),
list(label='PC4', values=~PC4),
list(label='PC5', values=~PC5)
),
text=~Taxa,
diagonal=list(visible=F),
marker = list(
color = ~color,
#colors = clist,
size = 5,
line = list(
width = 1,
color = 'rgb(1,1,1)'
)
)
)
fig11 <- fig11 %>%
layout(
title = "Scatterplot Matrix for PCA",
hovermode='closest',
autosize = F,
dragmode = 'select',
plot_bgcolor='rgba(1,1,1, 0)',
xaxis=axis,#list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=2, titlefont=list(size=13)),
yaxis=axis,#list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=2, titlefont=list(size=13)),
xaxis2=axis,
xaxis3=axis,
xaxis4=axis,
xaxis5=axis,
yaxis2=axis,
yaxis3=axis,
yaxis4=axis,
yaxis5=axis
)
fig11 <- fig11 %>% style(showupperhalf = F)
fig11
Plot fastStructure
#setwd("BOL21afZAXA00/")
#Make a list of paths and a list of .meanQ file names
ffiles <- list.files("fastStructure",full.names=T, pattern = "*.meanQ")
flist <- readQ(files=ffiles)
#imports the order of taxa saved in the phylo tree steps.
tree.order.names <- read.delim("order.csv",header = T,)$x
#sorts the taxa the same way for pop structure plots.
tree.order.inds <- match(tree.order.names,taxalabels)
flist <- sortQ(flist, "k")
run.names <- attr(flist,which = "names")
if(length(unique(sapply(flist,nrow)))==1) flist <- lapply(flist,"rownames<-",taxalabels)
indfun <- function(x) flist[[x]][tree.order.names,]
flist <- lapply(run.names,indfun)
attr(flist, which = "names") <- run.names
#lab1 <- read.delim("TreeGroups000.txt", header = FALSE)
##Make the Plots!
text.size <- 5
kl <- length(flist)
kl1 <- round(kl/2)
kind1 <- c(1:kl1)
kind2 <- c((kl1+1):kl)
IchooseK <- c(8:12)
#uncomment this to pick your own values!
kind1 <- IchooseK
fplotfilename1 <- "fplot1"
fplotfilename2 <- "fplot2"
fastStructurePlot1 <- plotQ(alignK(flist[kind1], type = "across"), showlegend = T, imgoutput="join",splab = as.character(kind1+1), showindlab=T, useindlab=T, returnplot=T, exportplot=T, basesize = text.size, imgtype="png", exportpath = "images/", height=2,width=44, sppos="left",dpi=600,outputfilename=fplotfilename1)
## Drawing plot ...
## images//fplot1.png exported.
#fastStructurePlot2 <- plotQ(alignK(flist[kind2], type = "across"), showlegend = T, imgoutput="join",splab = as.character(kind2+1), showindlab=T, useindlab=T, returnplot=T, exportplot=T, basesize = text.size, imgtype="png", exportpath = "images/", height=2,width=44, sppos="left",dpi=600,outputfilename=fplotfilename2)
#grid.arrange(fastStructurePlot1$plot[[1]])
#grid.arrange(fastStructurePlot2$plot[[1]])
#){width=65%}
fm1 <- image_read(paste("images/",fplotfilename1,".png",sep=""))
image_rotate(fm1, 90) %>% image_write(paste("../../",fplotfilename1,"rot90",".png",sep=""))
#fm2 <- image_read(paste("images/",fplotfilename2,".png",sep=""))
#image_rotate(fm1, 90) %>% image_write(paste("../../",fplotfilename2,"rot90",".png",sep=""))
knitr::include_graphics(paste(fplotfilename1,"rot90",".png",sep=""))
